run("../subroutines/MatlabDataConstructionCPI.m")
load ../data/DataforMatlab/IP.csv
load ../data/DataforMatlab/torn_CPI.csv

load '../data/DataforMatlab/cpi_pvec_tmp';
load '../data/DataforMatlab/UK_Budget_CPI';

I = size(B_vec,1);N = size(B_vec,2);T = size(B_vec,3);pvec = reshape(pvec_tmp',[I,1,T]);
[U_vec] = CalMoneyMetric(I_vec, B_vec, pvec, 0);
IP = IP(:,end-16:end);
 

%% Percentile for plotting
% Reshape data
% Reshape data
Wfull = reshape(I_vec(:,:,:), [N,T]);
U = reshape(U_vec,[N,T]);
% Initialize variables
PC = zeros(9,T);
Wq = zeros(9,T);
Uq = zeros(9,T);

for p = 1:9
    PC(p,:) = sum(Wfull <= IP(p,:));
    PC(PC == 0) = 1;
    Wq(p,:) = Wfull(sub2ind(size(Wfull), PC(p,:), 1:T));
    Uq(p,:) = U(sub2ind(size(U), PC(p,:), 1:T));
end

%%Insert NAN into income vectors for which the money metric cannot be calculated.
[~, indmax] = max(U);
[~, indmin] = min(U);
W = NaN(size(Wfull));
Wmax = Wfull(sub2ind(size(Wfull), indmax, 1:T));
Wmin = Wfull(sub2ind(size(Wfull), indmin, 1:T));
for tau = 1:T
    W(indmin(tau):indmax(tau), tau) = Wfull(indmin(tau):indmax(tau), tau);
end

logUmax = max(log(U));
logUmin = min(log(U));


%RC = [Wmin;Wq;Wmax]./exp(torn_CPI') * 52;
OurxaxisUK = [Wmin(:,end);Wq(:,end);Wmax(:,end)]*52;
OuryaxisUK = [exp(logUmin(:,end)+log(52)) ;(Uq(:,end)*52);exp(logUmax(:,end)+log(52)) ];



WelfareRelevant = log( [min(W)./min(U);Wq./Uq;max(W)./max(U)]);
diffp = WelfareRelevant - torn_CPI';

load ../data/DataforMatlab/agg_mmcpi.mat
%%
Floguagg = makefunction(log(agg_I_vec(:,:,T)),log(agg_U_vec(:,:,T)));
logUdiss = Floguagg(log(I_vec(:,:,T)));
indU = indmin(end):indmax(end);

hhw = hhw_vec_L(:,:,end);
Div = I_vec./reshape(exp(torn_CPI'),[1,1,T]);


%%
WelfareRelevant = log( [min(W)./min(U);Wq./Uq;max(W)./max(U)]);
axoptions={'scaled ticks = false',...
           'y tick label style={/pgf/number format/.cd, fixed, fixed zerofill,precision=0, set thousands separator={}}',...
           'x tick label style={/pgf/number format/.cd,precision=1, set thousands separator={}}',...
           'legend style={font=\normalsize}'}; %\scriptsize

figure('name','Figure O3 (b)','Position',[100 100 700 400],'NumberTitle','off')
bar(diffp(2:end-1,end)*100);hold on

set(gca, 'XTick', 1:9, 'XTickLabels', {'10%','20%','30%','40%','50%','60%','70%','80%','90%','highest'})

xlabel('Income Percentile in 2017')
ylabel('Bias')
ylim([-5 10]);
%matlab2tikz('../fig/bias_85goods.tex','extraAxisOptions',axoptions)

figure('name','Figure O3 (c)','Position',[100 100 700 400],'NumberTitle','off')
bar(agg_diffp(2:end-1,end)*100);hold on
set(gca, 'XTick', 1:9, 'XTickLabels', {'10%','20%','30%','40%','50%','60%','70%','80%','90%','highest'})
xlabel('Income Percentile in 2017')
ylabel('Bias')
ylim([-5 10]);
%matlab2tikz('../fig/bias_17goods.tex','extraAxisOptions',axoptions)


%% FigureO3 (a)
axoptions={'scaled ticks = false',...
           'y tick label style={/pgf/number format/.cd, fixed, fixed zerofill,precision=0, set thousands separator={}}',...
           'x tick label style={/pgf/number format/.cd,precision=1, set thousands separator={}}',...
           'legend style={font=\normalsize}'}; %\scriptsize

figure('name','Figure O3 (a)','NumberTitle','off')
h = loglog(OurxaxisUK ,OuryaxisUK,'LineWidth', 2); hold on
g1 = gca;
xLimits = g1.XLim;

loglog(agg_OurxaxisUK ,agg_OuryaxisUK,'LineWidth', 2);


xlim([agg_OurxaxisUK(1), xLimits(2)]);

ax = ancestor(h, 'axes');
ax.XAxis.Exponent = 0;
ax.XAxis.TickLabelFormat = '%.0f';
xticklabels({'10,000','100,000'})
yticklabels({'','10,000','100,000'})
xlabel('Income in 2017')
ylabel('2001 Income that gives the same utility as in 2017')
legend('Money metric (85 categories)','Money metric (17 categories)','Location','northwest')

%matlab2tikz('../fig/Money_dis.tex','extraAxisOptions',axoptions);



% Function to interpolate vectors with NaN values
function F_X = makefunction(X,Y)
    indexToKeep = ~isnan(X) & ~isnan(Y);
    X = X(indexToKeep); Y = Y(indexToKeep);    
    [X, indexUN] = unique(X);
    Y = Y(indexUN);
    F_X = griddedInterpolant(X,Y);
end
  